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We consider the Kolmogorov-Sinai entropy for dilute gases of A^ hard disks or spheres. This can 
be expanded in density as /iks cx nN\[ana'^ + B + 0{na'^) + 0{1/N)\, with a the diameter of the 
sphere or disk, n the density, and d the dimensionality of the system. We estimate the constant B 
by solving a linear differential equation for the approximate distribution of eigenvalues of the inverse 
radius of curvature tensor. We compare the resulting values of B both to previous estimates and 
to existing simulation results, finding very good agreement with the latter. Also, we compare the 
distribution of eigenvalues of the inverse radius of curvature tensor resulting from our calculations 
to new simulation results. For most of the spectrum the agreement between our calculations and 
the simulations again is very good. 
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I. INTRODUCTION 

It is generally believed that the approach to equilib- 
rium of a typical many-particle system, such as a gas or 
liquid will depend on its dynamical properties. Specif- 
ically, the more chaotic the system, the more rapidly 
its approach to (at least local) equilibrium will proceed. 
Furthermore, the apparent randomness of these systems, 
in spite of their fully deterministic microscopic behavior 
(at least for classical systems) has also been attributed 
to the chaotic nature of their dynamics. Discussions of 
this can be found e.g. in books by Dorfman [1| and Gas- 
pard y and review papers by Van Zon et al. [3, |j] . A 
very common and generally used measure of chaos is the 
Kolmogorov-Sinai entropy, which we will denote by /iks- 
In systems which are closed, the Kolmogorov-Sinai en- 
tropy ft-KS equals the sum of all positive Lyapunov ex- 
ponents, the average rates over very long times of diver- 
gence (or convergence) of infinitesimal perturbations. It 
describes the rate at which the system produces informa- 
tion about its phase-space trajectories, or equivalently 
about the distribution of density over phase space in 
some ensemble. In systems with escape, the Kolmogorov- 
Sinai entropy has also been connected to transport coef- 
ficients p-Q ■ In such systems it is no longer equal to the 
sum of all positive Lyapunov exponents. 

Chaotic properties such as the Lyapunov spectrum of 
systems of low |8l-[llj as well as high ;1.2j dimensionality, 
such as moving hard spheres or disks, have been stud- 
ied frequently. Extensive simulation work has been done 



on their Lyapunov spectra |13l - ll5| . and for low densi- 
ties analytic calculations have been done for the largest 
Lyapunov exp onent (a. [l6l - (l8| . the Kolmogorov-Sinai en- 
and for the smallest positive Lyapunov ex- 
20| . Analytic methods employing kinetic 



tropy [111. 



ponents [1 

theory have been applied to calculate chaotic properties. 
Agreement between analytic calculations and numerical 
results is generally good, but with respect to the KS- 
entropy there is one notorious exception, which is the 
central issue of the present paper. 

In this paper we consider a system consisting of N 
hard, spherical particles, of diameter a, at small number 
density n, in d dimensions {d = 2,3). We calculate the 
Kolmogorov-Sinai entropy in the low density approxima- 
tion, where it is expected to behave as [ll| 



hKs = NvA 



\n{na'^) +B + 0{na'^) 



-1^ 



(1) 



* A.S.deWijn@scienc e.ru.iiI] 
t ^H.vanBeijeren@uu.ril| 



The constant A has been calculated by Van Beijeren et 
al. in [lll|, but the results found there for B were unsat- 
isfactory. 

In this paper we present a more successful calculation 
of B, through the distribution of eigenvalues of the in- 
verse radius of curvature tensor. The calculation pre- 
sented here differs from that presented by De Wijn in 
Ref. [21|, in that it is far more elegant and less cumber- 
some and the agreement of the results with values found 
in simulations is better. On the other hand, the calcu- 
lation here is less systematic and it is not clear how to 
apply the results of this paper to calculating specific Lya- 
punov exponents, as can be done 22!| with the results of 
Ref. [m. 

The paper is organized as follows: In section llll we in- 
troduce Lyapunov exponents and review the properties 



of hard sphere dynamics in tangent space (the space in 
which the dynamics is described of infinitesimal devia- 
tions between nearby trajectories in phase space). In 
section IIIII we introduce the radius of curvature tensor 
and its inverse, relate the KS-entropy to the time aver- 
age of the trace of the inverse radius of curvature tensor 
and investigate the dynamics of these tensors both dur- 
ing free flight and at collisions. In llVl we present two ap- 
proximate calculations of the average distribution of the 
eigenvalues of the inverse radius of curvature tensor. In 
section FVl we compare the results of these calculations to 
those of numerical simulations and we also compare the 
resulting value for the coefficient B in Eq. ([T|) to those 
obtained in simulations and in previous calculations. Fi- 
nally, in section FVTl we present our conclusions. 




FIG. 1. Two particles at a collision in relative coordinates. 
The collision normal a is the unit vector pointing from the 
center of one particle to the center of the other. 



II. LYAPUNOV EXPONENTS AND DYNAMICS 
OF HARD SPHERES IN TANGENT SPACE 

This section is an abbreviated version of similar sec- 
tions in Refs. [20, [2l[. It appears here to make this 
paper more self-contained. For more details the reader 
may also consult Ref. [IJ]. Consider a system with an 
7\A-dimensional phase space F. At time i = the sys- 
tem is at an initial point 70 in this space. It evolves 
with time, according to 7(70, t). If the initial conditions 
are perturbed infinitesimally, by Sjo, the system evolves 
along an infinitesimally different path 7 -I- ^7, which can 
be specified by 



(57(7o,i) = M^„(i)-(57o 
with the matrix M^g{t) defined by 



'70 



w = 



dla 



(2) 



(3) 



The Lyapunov exponents are the possible average rates 
of growth or shrinkage of such perturbations, i.e., 



\i = lim -ln|^j(t)| , 

i— 7-00 t 



(4) 



where ^i{t) is the i-th eigenvalue of \\l\^^{t). For ergodic 
systems, the Lyapunov exponents are expected to be the 
same for almost all initial conditions. For each exponent 
there is a corresponding eigenvector of M^^{t). 

For a classical system of hard spheres without internal 
degrees of freedom, the phase space and tangent space 
may be represented by the positions and velocities of all 
particles and their infinitesimal deviations. 



7j = (rj,Vj) , 
^7i = {6vi,6wi) , 



(5) 
(6) 



where i runs over all particles and 7,; and 5^i are the 
contributions of particle « to 7 and ^7. 

In the case of a purely Hamiltonian system, such as 
the one under consideration here, hard spheres with only 



the hard particle interaction, the dynamics of the sys- 
tem are completely invariant under time reversal. To- 
gether with Liouville's theorem, which states that phase 
space volumes are invariant under the fiow, this leads to 
the conjugate pairing rule |2j, [2J|, i.e. for every posi- 
tive Lyapunov exponent there is a negative exponent of 
equal absolute value. In systems which are time rever- 
sal invariant, but do not satisfy Liouville's theorem, the 
conditions for and the form of the conjugate pairing rule 
are somewhat different p7l |. 

The system under consideration here has only hard- 
core interactions. Consequently, the evolution in phase 
space consists of an alternating sequence of free flights 
and collisions. 

During free flights the particles do not interact and the 
positions change linearly with the velocities. The compo- 
nents of the tangent-space vector accordingly transform 
to 






1 (i-io)l 
1 






(7) 



in which 1 is the d x d identity matrix. 

At a collision between particles i and j momentum 
is exchanged between the colliding particles along the 
collision normal, o- = (r^ — rj)/a, as shown in Fig. [T] 
The other particles do not interact. For convenience we 
switch to relative and center of mass coordinates, Srij — 
Sri — SvjjSIlij ~ {Sri + Srj)/2^_Svij= Svi — Svj, and 
SV,j ^ {Svi + Svj)/2. We find ^,^ 

(5r', — Sra — 2S ■ Sr^ , 



<5R',- 



'■V 

SR 



(8) 

(9) 

(10) 

(11) 

(12) 

Q ^ [(O- ■ V,j) 1 + (TV,j] ■ [(6- ■ Vij) 1 - Vij&] 

a(6--Vy) 



5v^ = Svij - 2S • Sv,j - 2Q ■ (Sr 
in which S and Q are the d x d matrices 



(Tcr 



where w^ = v^ — Vj . Here the notation ab denotes the 

Note that Q 

into vectors 

The vector v^ is a right 

v^ ■ a left zero eigenvector. 

determine M-,n(i). 



standard tensor product of vectors a and b, 
transforms vectors that are orthogonal to v,j 
that are orthogonal to v^ , 
zero eigenvector of Q, and 
Equations ^ through ^ 



'70 



III. RADIUS OF CURVATURE TENSOR 

Of particular interest for the Kolmogorov- Sinai en- 
tropy is the radius of curvature tensor, as the former 
is equal to the time average of the trace of the latter's 
inverse. Let 5y and 5w represent the full dN dimensional 
position respectively velocity perturbations of all parti- 
cles. The inverse radius of curvature tensor, T{t) is de- 
fined d, 1^, 1^ |251J26| as the inverse of the radius of cur- 



vature matrix 27]. It satisfies 

5v(t) = r{t) • <5r(t) , (14) 

resulting from some, almost arbitrary initial 7"(0) ^28j . 

A. KS-entropy and inverse ROC-tensor 

To establish the relationship between the KS entropy 
and the trace of the inverse radius of curvature tensor we 
consider the time evolution over long times of the pro- 
jection onto (5r space of the infinitesimal volume evolving 
from an initial volume in tangent space spanned by (5r(0) 
and (5v(0). Since the projections of all Lyapunov eigen- 
vectors with positive exponents onto 5y space are linearly 
independent, the size of this projected volume will grow 
roughly as the exponent of the sum of the positive Lya- 
punov exponents. A more precise statement is 



1 vol(5r(t) ■r-^ 
In — , . .„. = Z^ -^^ 

Ai>0 



lim 
t->oo t volor(O) 



(15) 



Here vol(5r(t) denotes the volume of the projection onto 
(5r space of the evolving infinitesimal volume in tangent 
space. Additionally, we have the identity 



5v(t) 



d6Y{t) 

dt 



(16) 



which holds, except at the instants of collisions t^, when 
(5r is refiected in a volume-preserving unitary transfor- 
mation Ui . Together with Eq. P^ , one obtains 

(t,/At-l 
JI^Yj lim W (l + r(7iAi)At)5r(0) 
i n=*i_i/At 

ti/At-l 

= ]JdetW,;lim ]J dct(l + r(ri Ai)At) , 



=ti_i/At 



where the first product is over the sequence of all colli- 
sions and we employ the convention YYi=i '^* = ^n • • • oi- 
This leads to 



,. 1, vol(5r(t) , _, 
lim -log ,.)1 ={^^T). 
t-s-oo t volor(O) 



(18) 



By comparing Eqs. (|T5t and ([T8| one obtains the desired 
identity ^ 



iKS 



(TrT) 



(19) 



Within the framework of this article, this identity is very 
central. It relates the trace of the inverse radius of cur- 
vature tensor directly to the Kolmogorov-Sinai entropy, 
which we wish to calculate. 

The rest of this paper is therefore dedicated to the de- 
scription of the eigenvalues of the inverse radius of curva- 
ture tensor. We will consider their dynamics and derive 
approximate equations for the time evolution of their dis- 
tribution. Using these we will calculate the changes in the 
distribution of the eigenvalues due to both free streaming 
and collisions. From this we obtain approximations for 
the stationary distribution of eigenvalues. 

This calculation is simplified appreciably by restric- 
tion to low densities. We may assume then that in each 
collision the elements of the precollisional inverse ROC 
tensor involving either of the colliding particles are small, 
due to their decrease during the free flight preceding the 
collision. This assumption is violated only in collisions 
where one of the colliding particles had collided shortly 
before. The fraction of all collisions where this is the case 
decreases linearly with density. 

In the next section two simple approximation schemes 
will be presented, based on the dynamics describing the 
changes in the distribution of eigenvalues of the inverse 
ROC tensor resulting from collisions. 



B. Dynamics of the inverse radius of curvature 
tensor 

At a given collision, between particles labelled i and 
j, let S and Q be the dN x dN dimensional matrices 
which perform the transformations of 25 and — 2Q on 
the relative components in tangent space of the colliding 
particles and act as zero on all other independent com- 
ponents of 5v or (5r, as is described in Eqs. ([SHTTj). Note 
that Q-{I — S), where I is the dN x dN identity matrix, 
has d — \ nonzero eigenvalues and is symmetric. One of 
^the non-zero eigenvalues of Q • (I — 5) is equal to 



^0 = 



2Vi 



(20) 



(17) 



with Wij the unit vector in the direction of Vjj. The 
corresponding eigenvector has components in the sub- 
spaces belonging to the colliding particles of e^ = — Gj = 



((Ty • vj,)vj ■ and ek — ioi k ^ i,j along the di- 



rections of the other particles. For d > 2 the other d — 2 
non-zero eigenvalues are given by 



(o 



2v,, 



a 



(21) 



with eigenvector with components e^ = — e^ normal to 
both Vij and & and again e^ = 0. The dynamics of the 
inverse radius of curvature tensor at a collision can be 



derived by expressing Sv' in terms of dr' , 



= [{i-s)-r-{i-sy 



(22) 



Q-{I-S)-^]5v' 



(23) 



With (X — iS) ^ = (I — 5) , we find for the inverse radius 
of curvature tensor after the collision, 

T = [I - S) -T ■ [I - S) + Q- [I - S) . (24) 

The dynamics during free flight follow from (fni) as 

T{t + dt) = [T{t)-^ + Idt]-^ (25) 

The eigenvectors of 7" do not change during a free flight, 
so its time evolution may be specified by the evolution of 
its eigenvalues ^. From Eq. (f25|) one finds these satisfy 



^'{t) = -m' 



with solution 



m 



^-1(0) +t 



(26) 



(27) 



Note that this may be simplified by considering the dy- 
namics of the radius of curvature tensor, T~^. For this 
operator, the drift velocity becomes a constant, 1, ir- 
respective of the eigenvalue and the choice of time unit. 
However, in this case Eq. (l23l) becomes more complicated. 



and left eigenvectors. The corresponding eigenvalues to 
linear order in n* are given by 



? 



[{i-s)-r- 



{I-S) + Q-{I-S)]-e, 



(29) 
(30) 



since the elements of T between these eigenvectors on 
average are of order n* compared to the elements of 
Q-{I-S). 

Under the approximation of Eq. (|30)) the vectors e^ 
only depend on the collision parameters Vij and it, as 
specified below Eqs. (PU)) and (HH). The remaining eigen- 
values of T' , to leading order in n*, can be identified as 
the eigenvalues of the projection VTV of the matrix T 
onto the dN — d+1 dimensional space orthogonal to the 
d — 1 eigenvectors e^ , as follows from standard perturba- 
tion theory. Hence it follows that they are interspersed 
between the precoUisional eigenvalues. This is worked 
out in Appendix A. 

These eigenvalues are distributed in roughly the same 
way as the eigenvalues of the full matrix T [30, [3l|. But, 
as the eigenvalues of VTP lie in between those of T, the 
distribution of these eigenvalues is slightly narrower than 
that of the eigenvalues of T- For more details on this see 
Appendix B. 

In the next section, we present two approximation 
schemes. In the first scheme, the narrowing will be ig- 
nored and the approximation will be made that the dis- 
tribution of eigenvalues of PTV is the same as that of T- 
The resulting equation for the distribution of eigenvalues 
can be solved analytically. Its solution is expressed in 
terms of the distribution /o(^o) of the non-zero eigenval- 
ues of Q • (I - 5). 

In the second scheme, a somewhat more refined ap- 
proximation is made, in which we assume that each eigen- 
vector that is changed significantly by the collision, but 
not created in it, can be written as a linear combination 
of exactly 2 precoUisional eigenvectors. Why this is an 
improvement will be argued in the discussion, where we 
will also briefiy discuss possibilities for further improve- 
ments. 



C. Eigenvalue distribution at low densities 



At low densities the mean free-flight time is given by 



1/9 



am) 



vo2n*Tr 2 



(28) 



with vq = {kBT/my^^ the thermal velocity and n* — 
na'^. The probability of a particle colliding within a 
free flight time of order a/vo (this is the typical order 
of ^"-^(0)) is of order n* and these events, to first ap- 
proximation, may be neglected. 

Let us denote a spanning set of the subspace spanned 
by the eigenvectors with non-zero eigenvalues of Q • (X — 
S) for a specific collision as ei through e^^.i. Because 
Q ■ {I — S) is symmetric, these vectors are both right 



IV. THE DISTRIBUTION OF EIGENVALUES 

OF THE INVERSE RADIUS OF CURVATURE 

TENSOR 

A. First approximation scheme 

Now that the dynamics of the eigenvalues of the inverse 
radius of curvature tensor are specified, we may write 
down approximate time-evolution equations for the dis- 
tribution of these eigenvalues, f{£,,t), which is normalized 
to unity. Let /o (1^) be the distribution of nonzero eigen- 
values of Q-{I-S), which follows from Eqs. (HO]) and ((2T|) 
and the distribution of the collision parameters v.y and 
(T. We equally normalize it to unity. The simplest ap- 
proximation for the rate of change of the distribution of 



eigenvalues of the inverse radius of curvature tensor is 
the one announced at the end of the previous section: af- 
ter a collision, the new distribution of eigenvalues is the 
same as the old one, except for the contributions from 
the nonzero eigenvalues oi Q ■ (T — S). Combining this 
with the rate of change resulting from free streaming one 
obtains 



dt 



fiU) 



i?{d-l) 
2d 



[m)-fiu)] 



d_ 



:efi^,t)] 



(31) 



where the single-particle collision frequency is given at 
low density by Eq. (1^5)) . The first term on the right-hand 
side is due to collisions. The first part of it is the gain. 
The second part is the loss. Its form here is based on 
our approximation that the shape of the distribution of 
remaining eigenvalues is not changed in a collision. The 
final term is due to the drift during free flight. 

For time going to inflnity, the distribution of eigenval- 
ues becomes stationary. In other words, the left hand 
side of Eq. (I5T]) becomes zero. Eq. ([51]) then may be 
rewritten in a more convenient form, as 

foiO = /stat(C) - cK'^/stat(e) + 2e/stat(0] , (32) 



d^ 



where 



2d 



(d-l)P ■ 
Solutions to this equation are of the form 



/stat (0 



1 

c'^o/o(6)-T2 exp 






(33) 



(34) 



Notice that for i^ small, as a consequence of Eq. (I^Ul) , this 
reduces to UAO = ^ exp (^) 

From Eq. ((T9)) it follows that the KS entropy may be 
obtained directly from the first moment of fstntiS,)- From 
Eq. dnH) we find that 



Nd 



Iks 



rfeo/o(Co)exp(^)r(0,^ 
c ./o \c4oJ V cfo 



(35) 



where r(, ) denotes the incomplete gamma function, de- 
fined by 



T{x,y) = / dte-^e' 
•'y 



Jo 



(36) 



At low densities the collision frequency is low, so that c 
is very large and the product exp (l/(c^o)) T (0, l/(c^o)), 
up to corrections of 0{n), is equal to In(c^o) ~ 7j where 
7 « 0.577216 is Euler's constant. The Kolmogorov-Sinai 
entropy then becomes 



^ND{d-l)(/ Co 
hKS^ 2 (1-7 



In 



2d 



(d-1) 



7 (37) 



We note that here the brackets, instead of a time aver- 
age denote an average over the probability distribution 
for the new eigenvalues at a collision, in this case /q. 
This can be expressed in terms of the joint probability 
distribution of the collision parameters as 



(giCo)) 



deo/o(eo)5(Co) 



/3m f d\ f , ^,., \-i 

^^^ V 2 j / «V'^Vj«'^^("Vy • '^)Kv« - Vj) ■ '^1 



</'M(Vi)'/'M(Vj 



-2v,, 



+ {d~2)g 



-2-Vij ■ or 



(38) 



where Eqs. (|20l) and ([?T|) have been substituted and 
</'m(v) is the Maxwell distribution. 



(/>m(v) 



/27rfcBr 

\ TO 



-d/2 



exp 



2fcBr 



(39) 



The function d(x) is the unit step function, which van- 
ishes for X < and equals unity for x > 0. In general, 
time averages of functions of ^ may be expressed as av- 
erages over fstat- 

In Sec. |V] the results from Eq. ([57]) will be discussed 
and compared with results from molecular dynamics sim- 
ulations. 



B. Second approximation scheme 

In the previous subsection, the distribution of eigen- 
values after a collision was assumed to be the same as 
the one before, except for the non-zero eigenvalues of 
Q • (X — 5) . In Appendix A it is shown that in reality 
these eigenvalues are determined by the equation 



E 



^-e. 



0, 



(40) 



at least in the case of d = 2, when Q ■ (T — 5) has a single 
non-zero eigenvector J32l | e, which can be expressed in 
terms of precoUisional eigenvectors ^i of T with eigen- 
values Ci, as e = ^iCitpi. From Eq. (|i(7| one sees that 
precisely one new eigenvalue ^ originates between each 
subsequent pair of precoUisional ones. 

We can divide the Ci into two categories, appreciable 
and almost vanishing. In a pragmatic way, this distinc- 
tion can be made by considering as appreciable the set 
of largest cf that sum to all but a small fraction of unity 
(for instance 0.01). The other cf , then, are almost van- 
ishing. Due to the locality of the interactions, one may 
argue that, in the limit of a large system, the number of 
appreciable Ci is small compared to N in most cases. 

For eigenvectors with very small values of Ci a new 
eigenvalue is found very close to an old one, typically to 



-"H^ 



-^ — ef- 



-^>f- 



D >o X 



FIG. 2. An impression of eigenvalues on an interval. 
The crosses denote precoUisional eigenvalues, the circles and 
squares postcoUisional ones. The bold and the regular crosses 
denote eigenvalues corresponding to eigenvectors with an ap- 
preciable respectively almost vanishing cf. The circles and 
square denote new eigenvalues. The eigenvalue indicated by 
the square is approximately the solution of Eq. (|40|) with 
only the eigenvalues with appreciable c^ , indicated with bold 
crosses, contributing. The circles almost coincide with prec- 
oUisional eigenvalues with almost vanishing cf . 



the left respectively the right of the old one if the nearest 
eigenvalue with appreciable Ci is to the left respectively 
the right of ipi . This may be interpreted in the following 
way: in Eq. (j40p one may ignore all eigenvectors with 
almost vanishing Ci, because their forms and eigenvalues 
remain essentially unchanged. For the remaining eigen- 
values one retains the property that new eigenvalues are 
interspersed between the old ones. This is illustrated in 
figure m 

The number of eigenvectors contributing appreciably 
to Eq. (HO)) varies from collision to collision, but it always 
is at least 2, because the two colliding particles cannot 
have collided before without intermediate collisions with 
other particles. The actual distribution of the number 
of contributing eigenvectors and the distribution of the 
values of the corresponding c^ are not easy to determine. 
In this subsection we make two simplifying assumptions: 
firstly that the number of contributing eigenvectors al- 
ways is just 2 and, secondly that their coefficients ci and 
C2 are distributed isotropically, irrespective of the eigen- 
values <^i and ^2- This means that these coefficients can 
be represented as ci = cos (f) and C2 = sin 0, with (f> dis- 
tributed uniformly on the unit circle. This assumption 
implies that the distribution of the corresponding original 
eigenvalues ^i and ^2 is the same as the (as yet unknown) 
overall distribution of eigenvalues of the inverse radius of 
curvature tensor. Obviously, these assumptions are at 
best approximately correct, but in our discussion we will 
make plausible why, on the basis of these assumptions, a 
better approximation can be obtained for the eigenvalue 
distribution than the one given in Eq. p4p . 

At a collision, the two eigenvalues ^1 and ^2 disappear 
and are replaced by a new eigenvalue, ^O: related to Q ■ 
(I — S), and the mixed eigenvalue 



iicl 



■6c?, 



(41) 



as follows from Eq. (PH)) . The eigenvector belonging to 
this mixed eigenvalue is a linear combination of the two 
old eigenvectors, orthogonalized to the non-zero eigen- 
vectors oi Q ■ {I — S). Under these approximations the 



collision term in Eq. ((3T|) is modified, leading to 



+§-^[ef{u)] 



(42) 



where fcoii[f]{S.,t) represents the distribution of the new 
mixed eigenvalue after the collision, as a functional of 
f{^,t), the distribution before the collision. 

Under the assumptions described above this distribu- 
tion can be written as 



fconUKtt) 



de'</(e',i)/(r,i)/i(eie,r)(43) 



Here /i(CI^'jC") is the distribution of the new eigenvalue 
^ between the eigenvalues £,' and £," . The function h 
assumes the form 

H^, \^', n = — ^ . (44) 



From this an equation similar to Eq. (j32p can be derived. 
One finds 

foiO — 2/stat(C) ~ /coiu/stat](0 

-cK'^/stat(C)+2C/stat(C)] , (45) 

This equation can easily be solved numerically and from 
its solution a second prediction of hus can be obtained. 
These results are discussed in the next section. 



V. RESULTS AND DISCUSSION 

In the previous sections we have developed two closely 
related analytical schemes that enable us to calculate ap- 
proximations for the stationary distribution of eigenval- 
ues of the inverse radius of curvature tensor. From this 
distribution we may obtain expressions for the leading or- 
der terms in the density expansion of the KS entropy of 
a gas of hard disks or spheres. In particular, the predic- 
tions resulting from Eq. p4p and numerical solutions of 
Eq. PSJ) may be compared to the results of Refs. [ll|, [21[ 
and results from molecular dynamics simulations. 

In Ref. [ll| Van Beijeren et al. proposed as approxi- 
mation for the KS-entropy 






Ni?{d-l) 



In Co + In 



(46) 



with Ti the free flight time of particle i since the previous 
collision. Putting this in the form of Eq. ([IJ leads to 
A ^ {d — l)/2 and, after numerical integration. 



B 



(0) 



0.209 



if rf = 2 



-0.583 if d == 3 



(47) 



From molecular dynamics simulations Posch and 
coworkers [111, [SJ] found the following results for the 
Kolmogorov- Sinai entropy at low densities: 



I, num. 
"KS 



(0.499 ±0.001)iVj>(- In na'' + 1.366 ±0.005) if d = 2 



(1.02 ± 0.02)ND (- hina'^ + 0.29 ± O.Ol) 



if d = 3 



(48) 



r 



Comparing the results of Eq. (|37| to those of Ref. [llj , 
Eq. (|46| . we find that A is the same, but for B one obtains 
corrections to Eq. (j47|) of the form 



AS = hi 



2d 



1 



7- 



hi 



vin + Tj) 



(49) 



Dorfman et ah aheady expected corrections of In 4 w 
1.386 for d = 2 and hi 3 « 1.098 for d = 3 34], corre- 
sponding to the first term in Eq. (HHl) . 

In Ref. (2l| , elements of the radius of curvature matrix 
were estimated by considering the stretching of the tan- 
gent phase space during a sequence of two collisions with 
free flights, it was estimated that 



B- 



dW 



1.47 ±0.11 
0.35 ±0.08 



if d- 
if d- 



(50) 



We have evaluated the averages in Eq. ([^^ by integrat- 
ing over the joint distribution of the collision parameters 
[see Eq. (l38l) ]. The values for the parameter B resulting 
from the first approximation scheme follow as 



5(1) = 



7 + In 2 



iluTT 



|7 + ln3- iluTT 



1.255 
: 0.160 



if d= 2 
if d = 3 

(51) 



These results are in reasonable agreement with the results 
from the molecular dynamics simulations [33,] . given in 

Eq. dm). 

More accurate results can be obtained from the second 
approximation scheme by numerically solving Eq. (|45p. 
The solution for /stat(C) for d = 2, n = 0.001 is displayed 
in Fig. |4l From this, one finds an additional correction 
to B of 



AB = 0.086 



(52) 



regardless of dimensionality. This leads to a final result 
for the constant B of 



5(2) 



1.341 if d = 2 
0.247 if d = 3 



(53) 



This is in good agreement with the results from the 
molecular dynamics simulations, and in particular also in 
better agreement than the results of Ref. 21], Eq. ((50)) . 



A. Comparing approximation schemes 

We now argue why the distribution of inverse radius 
of curvature tensor eigenvalues obtained from the two- 
eigenvalue approximation resulting into Eq. (I42[) can be 



expected to be better than the simpler approximation, 
Eq. ([5T|) , resulting from assuming the distribution of in- 
terspersed new eigenvalues to be the same as that of 
the precoUisional eigenvalues. Consider the first two 
moments of these distributions. In the simple approx- 
imation these are not changed from their precoUisional 
values. Therefore, as mentioned already, the spectrum 
does not exhibit any narrowing in collisions, as it should 
according to the arguments presented before, which are 
supported by the calculations presented in Appendix B. 
In the two-eigenvalue approximation, two eigenvalues ^i 
and ^2 are sampled independently from the stationary 
distribution and replaced by one interspersed eigenvalue 
with the value £,' = Cici + £,2c1, according to Eq. (|41l) . 
The other interspersed eigenvalues retain their precol- 
lisional values. The coefficients ci and C2 are sampled 
as ci = cos (j) and C2 = sin (f>, with (p distributed uni- 
formly on the unit circle. Hence the average value of 
^' is the same as that of the precoUisional eigenvalues, 
as should be the case (see Appendix B). From Eq. (|4T|) 
and the assumed distribution of the Ci one also easily 
finds the coUisional changes of the second moments. In 
Appendix B it is shown that, to leading order in 1/n 
the average of the second moment is reduced at a col- 
lision by a factor 1 — (1 — A)/n, with A a constant 
with a value between zero and unity. This constant, de- 
fined in the Appendix through J2^=i i^t (^« ^ iO) ) = 

A/n (J2^=i id ~ {Oy)i '^^^ be calculated in the two- 
eigenvalue approximation from the assumed distribution 
of ci and C2 as A2 = 3/4. Since in most cases the new 
eigenvector will be composed of more than 2 precoUi- 
sional eigenvectors, A2 will be an upper bound to the 
actual A. Hence the two-eigenvalue approximation does 
lead to a narrowing of the spectrum, but it underesti- 
mates its extent. This is especially true in the region 
where /(^) reaches its maximum, since there the eigen- 
functions of the ROC tensor tend to be carried by many 
particles, as one can see from Figs. [3] and ID 

The results may be improved, in principle, by con- 
sidering larger sets of eigenvectors for spanning the e^. 
However, to do this in a sensible way one would need 
the distribution of the Ci, preferably as a function of all 
^i. So far no theory has been developed for this and it 
seems no simple task to do so. One could of course study 
this distribution numerically, but that would bring one 
close already to a full numerical study of the eigenvalue 
spectrum of the inverse radius of curvature tensor. 

In Fig. [3] we plot numerical results for the average 
number of particles contributing to an eigenvalue as 
function of ^, in a system of 64 two-dimensional par- 
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FIG. 3. Eigenvalue composition, as obtained from molecular 
dynamics simulation of 64 two-dimensional particles at n* 
0.01 in a square box. The figure shows an estimate for the 
average number of particles which contribute significantly to 
an eigenvector of the inverse radius of curvature as a function 
of the corresponding eigenvalue ^. This estimate is given by 
the ratio of the averages of the second and fourth powers of 
the component of the eigenvector in the subspace of a particle. 
The estimate is somewhat indirect, but it is clear that it is 
exact in the special case that the weight is distributed equally 
over a set of particles. At this density, the collision frequency 
is approximately equal to 0.0354(fcBr/m)'^"/a. In order to 
obtain sufficient information both at high and low densities, 
non-linear binning was used. The plot only shows results for 
bins that contain more than 10 points. 



FIG. 4. The distribution of eigenvalues of the inverse ra- 
dius of curvature tensor for d = 2 as calculated in this paper 
compared to results from simulations for various densities. 
Dashed and dot-dashed lines show the results of the present 
calculations, while solid lines represent the values found in the 
simulations. The results for the second approximation scheme 
are only shown for n* = 0.001, for which we calculated them 
with great accuracy. Similar curves would be obtained at the 
other densities. The simulated systems consist of 64 parti- 
cles in square boxes with periodic boundary conditions. Runs 
were performed of 100000 collisions, and the eigenvalues of 
the radius of curvature were calculated and binned every 64 
collisions. In order to obtain sufficient information both at 
high and low densities, non-linear binning was used. The plot 
only shows results for bins that contain more than 10 points. 



tides in a square box with periodic boundary condi- 
tions. For ^ larger than the collision frequency most 
eigenvectors are carried by two particles, indicating that 
these are new eigenvectors e.;. For smaller ^ there is a 
rapid increase in the average number of particles car- 
rying an eigenvector, followed by a sharp drop below 
^ w 0.25i/ « 0.9 X 10-2(fcBT/m)^/Va- The eigenvectors 
corresponding to smaller eigenvalues have drifted for a 
long time, during which in most cases the particles con- 
tributing to them have collided many times. These eigen- 
vectors are therefore typically carried by more particles. 
But, remarkably, for very low eigenvalues the number of 
particles carrying the eigenvector becomes very close to 
1. This is related to the existence of particles that have 
not collided for several mean free times. As a result of 
subsequent projections normal to new eigenvectors the 
weight of such a particle in the remaining eigenvector 
can increase from the original value 1/2 to values close 
to unity. The contribution of these eigenvectors to the 
KS-entropy is very small. 

We cannot directly translate the data of Fig. [3] to an 
estimate of the number of eigenvectors contributing sig- 
nificantly to a newly generated eigenvector e^. It is clear 
though that there is a strong correlation. Since the new 
Ei are always carried by just the two colliding particles. 



the components along them of eigenvalues carried by sev- 
eral particles necessarily have to be small. Hence, many 
of these are required to reconstruct any given e^ . 

Another, more technical approach was based on a cal- 
culation of the distribution of elements of the radius of 
curvature [2l|, [S^ • The results in the present paper are 
more accurate than the results of that calculation, and 
were obtained in a more elegant way. On the other hand, 
it is not directly clear to us how to further improve the 
accuracy of the calculation presented here, nor if the dis- 
tribution of eigenvalues of the radius of curvature could 
be used to calculate specific Lyapunov exponents of the 
system, as can be done with the distribution of the ele- 
ments of the radius of curvature 1231 ■ 



B. Distribution of eigenvalues from simulation 

In our present calculation of the Kolmogorov-Sinai en- 
tropy of a dilute hard-sphere gas the central quantity to 
be computed is the stationary distribution of eigenvalues 
of the inverse radius of curvature tensor. It is interest- 
ing to compare the calculated distribution to results from 
computer simulations. We have performed MD simula- 
tions for a system of hard disks, in which we calculated 



the radius of curvature tensor from the numerical values 
of 6r and 6v by making use of Eq. ([T^ and diagonalized 
it at regular time intervals. The results of these simula- 
tions are displayed in Fig. 21 along with the theoretical 
predictions. For a large range of eigenvalues, the calcu- 
lations follow the simulations closely, including at very 
high ^. At the lowest values of the range studied, there 
appear some differences. 

It can be seen from Fig. U that the small eigenvalues 
have a more peaked, hence narrower distribution than 
was found from the calculations. This is due to the fact 
that only linear combinations of two eigenvectors were 
considered for e^. In fact, as can be seen from Fig. [3] 
eigenvectors near the peak are generally carried by more 
particles and their contribution to e^ , which is carried by 
only two particles therefore will have a small coefficient 
Ci. As argued above, the larger the number of parti- 
cles carrying an eigenvector, the smaller the value of A 
and the stronger the narrowing. And Eq. (I40[) reveals 
that for given ^ the dynamics are dominated by nearby 
eigenvalues and their q. As can be seen from Fig. 21 
the distribution calculated using linear combinations of 
two eigenvectors not only predicts the KS-entropy better 
than the approximation based on no change in the spec- 
trum of interspersed eigenvalues, but it also follows the 
simulation results more closely for intermediate eigenval- 
ues. 



eigenvectors. This is highly nontrivial, however. 

It should be noted that, though the specific dynam- 
ics of the inverse radius of curvature tensor are differ- 
ent for other high-dimensional systems, such as the high- 
dimensional Lorentz gas 12] (which has uniformly con- 
vex scatterers), their overall behavior is generic for all 
systems consisting of many particles. 

An attractive approach seems describing the dynamics 
of the eigenvalues of the inverse radius of curvature ten- 
sor by means of a Fokker-Planck equation. In order to 
do so one needs expressions for the local drift and diffu- 
sion of the eigenvalues due to collisions. For this again 
more knowledge is needed of the way new eigenvectors are 
composed of old ones. Furthermore, the Fokker-Planck 
equation is a good approximation for systems where the 
dynamics consist of small jumps, while in the present case 
also large jumps happen. However, these mostly will oc- 
cur for large values of ^, where the dynamics is dominated 
by the drift. Therefore the Fokker-Planck equation may 
still be a good approximation. 

Finally, the expressions derived here for the dynamics 
of the inverse radius of curvature tensor, and the equa- 
tions for the distribution of its eigenvalues, Eqs. ([3T|) 
and (UH), can also be used for systems in a stationary 
non-equilibrium state. In such a state, the distributions 
of velocities and collision parameters are different and 
one has to take this into account when calculating the 
averages or the source terms in Eqs. (PT|) and (|^. 



VI. CONCLUSIONS 

In this paper we have calculated the Kolmogorov-Sinai 
entropy of systems consisting of hard disks or spheres 
from the stationary distribution of the eigenvalues of the 
inverse radius of curvature tensor. The dynamics of these 
eigenvalues consist of free streaming and collisional ef- 
fects. The latter are a combination of the generation 
of new eigenvalues, with a well-defined distribution, and 
a slight narrowing of the spectrum of remaining eigen- 
values with respect to the precollisional spectrum. Al- 
ready a simple approximation, which ignores this nar- 
rowing and assumes the spectrum to remain unchanged 
at collisions on average, reproduces the numerically ob- 
served spectrum quite well, with a fairly accurate predic- 
tion of the KS entropy. A slightly more refined approx- 
imation, assuming the new eigenvectors consist of just 
two precollisional ones, both sampled randomly from the 
full distribution, does predict a narrowing of the spec- 
trum at collisions, though it underestimates its extent. 
This approximation gives quite accurate results for the 
KS entropy and it reproduces the eigenvalue spectrum 
of the inverse radius of curvature tensor better than the 
simplest approximation. The remaining underestimation 
of the narrowing is strongest and most clearly visible at 
small eigenvalues. 

In order to improve on the estimates developed here 
one needs more knowledge on the decomposition of the 
new eigenvectors created at a collision into precollisional 
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Appendix A: Eigenvalues of a submatrix 

Let Vi through Vd-i be the projection operators which 
project out the d — 1 orthogonal unit vectors ei through 
ed-i- We are interested in the (distribution of) nonzero 
eigenvalues of submatrix VTV of a symmetric matrix T, 
with V ~Vi- ■ ■ .-Vd^i- We shall determine these by con- 
sidering the nonzero eigenvalues ^'^ of T' = VxTVx, and 
their corresponding eigenvectors rji. The eigenvalues of 
VTV can now be determined by applying this procedure 
d— 1 times projecting subsequently normal to each of the 
d — 1 non- vanishing eigenvalues of (I — <S) • Q. Thus one 
obtains matrices T with a decreasing number of nonzero 
eigenvalues. 

Let "01 through i/'z be the z normalized eigenvectors 
of a matrix 7" with corresponding eigenvalues t,i. Let us 
write the unit vector that is to be projected out, e^;, and 



an eigenvector 77 of T' with eigenvalue £,' in terms of the 
eigenvectors of T, 



with 



i 

i 



(Al) 
(A2) 

(A3) 



As 7] has no component along e^, we may write 

ev^eJ2l^^^^ (A4) 

i 

i 
i 

with ^ a constant such that 



By taking the inner product of Eq. (|A6I) with a given '(/'i 
one finds that, 



A = -77 



MC* 



e-e. 



(A8) 



By substituting Eq. (lASp into Eq. (jA7p , and dividing by 
II and ^', we find that 



E 



6-e' 



. 



(A9) 
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From this equation it follows directly that between each 
subsequent pair ^i and Ci+i there must be precisely one 
solution for ^'. 



Appendix B: Narrowing of eigenvalue spectrum 

In order to investigate the narrowing of the spectrum 
of eigenvalues as result of a collision we rewrite Eq. (jA9l) , 
multiplying it by — Jli(^' ~ d) ^^id find 

e'-' E(i - ^')e^^'""' + E(i ~^"- '^'M^r-' 

i i<j 

+ •■• = 0, (Bl) 

where we have made use of Eq. (jA3p . By comparing the 
coefficient of (C')"~^ to the coefficient in the eigenvalue 
equation for f, one immediately obtains 



E^^-E(i-^')^» 

1=1 i=l 



(B2) 



In section IIIII we have found that the the eigenvectors 
Ez consist of equal and opposite components on two 
arbitrarily determined (colliding) particles along a unit 
vector that is distributed isotropically in d-dimensional 
space. From this it follows immediately that the average 
of c^ over many collisions has to be equal to 1/n. If in 
Eq. (|B2p we replace cf by this average we find that the 
mean value of the interspersed eigenvalues on average is 
the same as that of the precoUisional ones. Note that 
Eqs. (I5T]) and (H^ both satisfy this property. 

Similarly, from the coefficient of (C')"~^ '^^^ obtains 
the identity 



E 



?Z?7 



E(i-^? 



c?)60 



(B3) 



Combining this with Eq. (|B2[) one finds the identity 



n-l 



E m - (Of) = 4 ( E(e^ - (0) ) + ( 1 - ! - i ) E((e^ - (0)^> + T.i'^tid - (Of) 



(B4) 



Here the brackets indicate an average over many subse- 
quent collisions. We used the identity 



and we introduced the symbol (^) defined by 

1 " 



(B5) 



(B6) 



I 

To leading order in 1/n the terms proportional to 1/n^ 
in Eq. (jB4[) may be ignored. We introduce the constant 
A defined through 



A 



E(^'(^^-(0)') = -E((^^-(0) 



(B7) 



i=l 



Note that A is smaller than 1, since cj < cf and < cf >= 
1/n. To leading order in 1/n combination of Eqs. (jB4[) 
and (|B7p leads to the reduction factor mentioned in sec- 
tion lV Al Obviously, the smaller A, the stronger the nar- 
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rowing. This should apply also locally, implying stronger 



narrowing in regions where the eigenvectors tend to be 
carried by more particles. 
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